/*
时光荏苒
 * 引用命名空间：
 * using System;
 * using System.Numerics;
 * using System.Threading.Tasks;
 * 
 * .Net版本：.NET Core 2.1
 * 
 * 作者：屈永虎
 * 电话：18372628983
 * QQ：1226950821
 */
using System;
using System.Numerics;
using System.Threading.Tasks;

public class Shiguangrenrang_Programme
{
    /// <summary>
    /// 矩阵乘法
    /// </summary>
    /// <param name="left">第一个矩阵</param>
    /// <param name="right">第二个矩阵</param>
    /// <returns>结果矩阵</returns>
    public static double[,] MatrixMult(double[,] left, double[,] right)
    {
        if (left.GetLength(1) != right.GetLength(0))
            throw new Exception("第一个矩阵的列数和第二个矩阵的行数不同");
        int crNum = left.GetLength(1);
        int offset = Vector<double>.Count;
        int verAryLength = (left.GetLength(1) - 1) / offset + 1;
        int rowCount = left.GetLength(0);
        int colCount = right.GetLength(1);
        double[][] resAry = new double[rowCount][];
        Vector<double>[][] rowAry = new Vector<double>[rowCount][];
        Vector<double>[][] colAry = new Vector<double>[colCount][];

        ParallelOptions DefaultOptions = new ParallelOptions() { MaxDegreeOfParallelism = Environment.ProcessorCount * 2 };
        Func<double[,], int, int, double> getRowValue = new Func<double[,], int, int, double>((ary, idx0, idx1) => { return ary[idx0, idx1]; });
        Func<double[,], int, int, double> getColValue = new Func<double[,], int, int, double>((ary, idx0, idx1) => { return ary[idx1, idx0]; });
        Action<Vector<double>[][], int, double[,]> initMatAry = new Action<Vector<double>[][], int, double[,]>((matAry, dim, ary) =>
        {
            Func<double[,], int, int, double> getValue = (dim == 0) ? getRowValue : getColValue;
            int dimCount = matAry.Length;
            int count = verAryLength - 1;
            Parallel.For(0, dimCount, DefaultOptions, index =>
            {
                Vector<double>[] verAry = new Vector<double>[verAryLength];
                double[] temp = new double[offset];
                for (int i = 0; i < count; i++)
                {
                    for (int j = 0; j < offset; j++)
                    {
                        temp[j] = getValue(ary, index, j + i * offset);
                    }
                    verAry[i] = new Vector<double>(temp, 0);
                }
                Array.Clear(temp, 0, offset);
                for (int i = count * offset; i < crNum; i++)
                {
                    temp[i % offset] = getValue(ary, index, i);
                }
                verAry[count] = new Vector<double>(temp, 0);
                matAry[index] = verAry;
            });
        });
        initMatAry(rowAry, 0, left);
        initMatAry(colAry, 1, right);

        Parallel.For(0, rowCount, DefaultOptions, rowIndex =>
        {
            double[] row = new double[colCount];
            for (int colIndex = 0; colIndex < colCount; colIndex++)
            {
                Vector<double> vres = Vector<double>.Zero;
                Vector<double>[] rowVerAry = rowAry[rowIndex];
                Vector<double>[] colVerAry = colAry[colIndex];
                double dres = 0;
                for (int i = 0; i < verAryLength; i++)
                {
                    vres += (rowVerAry[i] * colVerAry[i]);
                }
                for (int j = 0; j < offset; j++)
                {
                    dres += vres[j];
                }
                row[colIndex] = dres;
            }
            resAry[rowIndex] = row;
        });

        double[,] resultMat = new double[rowCount, colCount];
        Parallel.For(0, rowCount, DefaultOptions, rowIndex =>
        {
            double[] row = resAry[rowIndex];
            for (int colIndex = 0; colIndex < colCount; colIndex++)
            {
                resultMat[rowIndex, colIndex] = row[colIndex];
            }
        });

        return resultMat;
    }
}
